A combined R-matrix eigenstate basis set and finite-differences propagation method 
for the time-dependent Schrodinger equation: the one-electron case 
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In this work we present the theoretical framework for the solution of the time-dependent 
Schrodinger equation (TDSE) of atomic and molecular systems under strong electromagnetic fields 
with the configuration space of the electron's coordinates separated over two regions, that is regions 
I and II. In region / the solution of the TDSE is obtained by an R-matrix basis set representa- 
tion of the time- dependent wavefunction. In region II a grid representation of the wavefunction 
is considered and propagation in space and time is obtained through the finite-differences method. 
It appears this is the first time a combination of basis set and grid methods has been put forward 
for tackling multi-region time-dependent problems. In both regions, a high-order explicit scheme 
is employed for the time propagation. While, in a purely hydrogenic system no approximation is 
involved due to this separation, in multi-electron systems the validity and the usefulness of the 
present method relies on the basic assumption of R-matrix theory, namely that beyond a certain 
distance (encompassing region /) a single ejected electron is distinguishable from the other electrons 
of the multi-electron system and evolves there (region II) effectively as a one-electron system. The 
method is developed in detail for single active electron systems and applied to the exemplar case of 
the hydrogen atom in an intense laser field. 



I. INTRODUCTION 

Exploration of the fundamental processes that occur 
when atomic and molecular systems are subject to ex- 
treme conditions is currently a major research area. 
Experimentally, such processes are realized by strong 
and/or short intense laser pulses radiating at infrared 
wavelengths 0, Q and have recently been utilised at a 
more practical level for reconstruction of nuclear proba- 
bility distributions, visualisation of molecular orbitals, 
alignment of molecules as well as production of high- 
order harmonics which in turn are used for the gen- 
eration of ultra short fields at the attosecond scale 

Theoretically, it is a huge task to treat the exact time- 
dependent (TD) response of a multi-electron system sub- 
ject to a strong electromagnetic (EM) field by ab initio 
methods. In response to extensive experimental achieve- 
ments using high-intensity Ti:Sapphire laser sources in 
the long wavelength regime, many theoretical studies em- 
ployed the strong-field approximation where the influ- 
ence of the Coulomb potential on the ejected electron 
wave function is neglected in favour of the external field. 
A more sophisticated approach that adopts the single- 
active-electron (SAE) approximation was also applied to 
the atomic case [9j. SAE models where one reduces the 
dimensionality of the multi-electron problem by freez- 
ing the most tightly bound electrons have proven to be 
very useful in cases where multiple electronic excitations 
are insignificant, and the SAE approximation is probably 
the most widely used approach when studying phenom- 
ena such as single ionization, above-threshold ionization 
(ATI) and high-harmonic generation (HHG). 

For systems of only two electrons, such as the negative 



hydrogen ion, helium, molecular hydrogen, direct, ab- 
initio, solutions of the time-dependent Schrodinger equa- 
tion (TDSE) have appeared in the early nineties (for a re- 
view see ref. [HI). Since then, the computational power 
has increased steadily and as a result these methods have 
reached a high level of accuracy, efficiency and reliabil- 
ity, tackling successfully the very demanding theoretical 
problem, of single and double ionization of helium at 390 
and/or 780 nm [HEl. 

Recently, the construction of FEL sources which de- 
liver brilliant radiation in the soft- and (in the immediate 
future) hard X-ray regime have initiated new challenges 
in the field of atomic and molecular physics [H, [Til ]. 
However, in contrast to what occurs with conventional 
laser sources, more than a single electron at a time re- 
sponds to short wavelength FEL light and X-ray FEL 
light will interact preferentially with the inner-most elec- 
trons, residing closer to the system's core, rather than 
with the valence ones. An immediate consequence of 
the above property is that theories such as the SAE and 
models not taking into account interelectronic interac- 
tions at a sufficient level are inadequate to describe the 
processes involved. Moreover, high-order harmonic gen- 
eration (HOHG) techniques are nowadays able to create 
pulses of subfemtosecond duration. Given that relaxation 
processes, such as Auger transitions, of the bound elec- 
trons are of the order of a femtosecond or less it can be 
concluded that the short time-variation of the EM field 
requires approaches where multi-electron dynamics can 
be reliably described. 

Given our intention to study multi-electron systems 
under intense EM ultrashort fields, there is consider- 
able importance in the development of computationally 
tractable methods able to treat multi-electron systems 
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with the least approximations possible. Such approaches 
have been developed in atomic and molecular physics 
studies, and include variants of time-dependent Hartree- 
Fock (TDHF) Q. Though a vast number of theoretical 
efforts in the spirit of TDHF [H, (M [H El El have ap- 
peared, even some extensions to include correlation be- 
tween the electrons, the question of how much and un- 
der what conditions correlation beyond the Hartree-Fock 
model is important still remains unanswered. The under- 
lying reason is the difficulties introduced by the nonlinear 
nature of the TDHF equations in combination with the 
fact that the single-configuration ansatz and the excita- 
tion process induced by the EM field are inconsistent. 
Improvements of the restricted Hartree-Fock ansatz and 
inclusion of exchange effects appear to be possible so- 
lutions to overcome such problems, although the appli- 
cations so far are only in one-dimensional (ID) models 

EE [M E3, m, m . 

An alternative ab-initio approach capable of treating 
multi-electron systems is R-matrix theory, with the ba- 
sic formulation appeared first in the context of nuclear 
theory, and later on applied in the field of atomic physics 
( [IE, Ea H3] ) ■ Traditionally, R-matrix theory is a theory 
where time is not involved in the study of the collision or 
photoionization processes. Variants of R-matrix theories 
and computational codes have been applied to an impres- 
sive number of systems, over the last 40 years (H)]. With 
the advent of strong and/or short laser pulse technology 
an early application of R-matrix theory to multiphoton 
processes appeared in the form of a Floquet expansion of 
the driven time-dependent wavefunction [29]. Although 
able to treat the field non-perturbatively, the R-matrix 
Floquet approach cannot be considered as a fully TDSE 
solution methodology since it is only suited to laser pulses 
containing many cycles. 

Similarly the appearance of high power sources at the 
short wavelength regime has led a number of theoretical 
groups to develop TDSE approaches based on R-matrix 
theory ( [30L l3ll l32| ) . with the first work to this end ap- 
pearing some years ago [Hj]. The basic assumption of 
R-matrix theory is very well suited to the physical sit- 
uation of the photoionization process involved in light- 
matter interaction. Under strong radiation any system 
will ionize either multiply or singly. In the regime of 
single ionization the ejected electron, after some time, 
depending on its distance from the core, can be safely 
identified as distinguishable from the other electrons. In 
R-matrix theory this is taken into account through the 
division of configuration space into two regions where, 
in the inner region (region /), all interelectronic and ex- 
change effects between all the electrons are treated, while 
in the outer region (region II) the ejected electron evolves 
effectively as a one-electron system under the influence 
of the residual core and the potential due to the remain- 
ing electrons. Thus in the outer region no matter what 
particular process has taken place the system wavefunc- 
tion consists entirely of that of the wavefunction of the 
ejected electron. 



The purpose of this work is two-fold. The first is 
to pursue development of a method which meets the 
above requirements for more complex systems than one- 
and two- electron systems and where atomic structure 
plays an important role in the processes. For this, a 
method based on R-matrix basis eigenstates appears to 
be tractable due to its success in describing such complex 
systems. Second, and equally important is the issue of 
efficiency and accuracy. It is inevitable that the demands 
of the calculations will make the study of such problems 
computationally very demanding. Finite-differences with 
high-order explicit time propagators [34j although diffi- 
cult to use throughout configuration space in a direct 
extension to multi-electron systems, have proven to be 
very efficient and accurate in solving the TDSE for one- 
and two-electron systems. In fact the HELIUM code [34| 
using such methods to solve the TDSE fully for a two- 
electron atom exposed to intense laser fields is able to 
run with high efficiency in both computation and com- 
munication over many thousands of cores on the largest 
supercomputers presently available. This established ef- 
ficiency makes their implementation for the outer region 
in our present approach a very reasonable one. In region 
/, an R-matrix basis set is used to propagate the multi- 
electron wavefunction while in region II amounting effec- 
tively to a one-electron problem, a finite-difference high- 
order propagation algorithm is used. Since to the best 
of our knowledge no such attempt has appeared, namely 
the propagation of the TDSE in a combined basis and 
grid representation of the TD wavefunction, we consider 
it essential to set out carefully in detail the basics of 
the method, free from complications arising from multi- 
electron considerations. Thus, we provide below the de- 
tails of the method and its usefulness for one-electron sys- 
tems and present results for the hydrogen system where 
accurate ab-initio methods, to compare with, are avail- 
able to us. 

The paper is organized as follows. In Sec. [IT] we give 
an overview of the basic ideas and principles. Section Unl 
is the key section of this paper and there we set out in de- 
tail the theoretical formulation for a one electron system. 
In Sec. IIVI we apply the method to the hydrogen atom 
in an intense laser field which serves as an exemplar. We 
have relegated to appendices some of the more technical 
details. Finally we set out some conclusions and perspec- 
tives with regard to the new method in Sec. [V] Atomic 
units are used (m = % = \e\ = ao = 1) throughout. 



II. OVERVIEW OF THE BASIC IDEAS AND 
PRINCIPLES 

As mentioned briefly in the introduction, the basic as- 
sumption of R-matrix theory for the outer-region wave- 
function allows the derivation of a TDSE (in the outer 
region), where only one electron is involved reducing 
the dimensionality of the problem there to its minimum, 
namely to at most three, thus simplifying the compu- 
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tational problem considerably. To put this in a more 
quantitative fashion, let us recall the (N + 1)— electron 
wavefunction beyond a certain distance, say b (taken as 
the inner boundary of the outer region II) (28j : 



Finite differences / R-matrix basi 



ijj(f N ,r;t) = y^$ 7 (fjv;f,crAr +1 )i/ 7 (r,i) r > b, 
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(1) 

with Fat = (n, r2, .., rjv), Ti < b,i = 1,2, ..,N and r — 
rjy+i. The < & 7 (rjv; f, ctat+i) are channel functions formed 
by coupling the target states of the residual atomic sys- 
tem ^ 7 (?jv), described by the Hamiltonian iJ/v(?/v) an d 
the angular and spin quantum numbers of the ejected 
electron. The radial motion of the ejected electron (in 
the 7-channel) is described by the radial channel func- 
tions / 7 (r, t) . The absence of the antisymmetrization op- 
erator is essential in the above expansion since it relies 
on the ejected electron and the remaining iV— electrons 
occupying different portions of configuration space, thus 
making the ejected electron distinguishable from the oth- 
ers. Let us now consider the TDSE of the above system, 
in an external time-dependent radiation field. By writ- 
ing the Hamiltonian for the field- free (N + 1)— electron 
system as -ff (f N , r) = + H N (r N ) + V(tn, r) we 

end up with the following form for the TDSE: 

d 

(2) 

with Z?(?7v,r,t) denoting the interaction operator be- 
tween the system and the external field, in the dipole 
approximation. Projection of the known channel states 
<£> 7 onto the TDSE and integration over Fjv and f, vn+i 
results in the following set of coupled partial differential 
equation for the radial motion in channels 7, 



d 



/ 7 (r, t) = krtffyir, t) + ^ D 7>r (r, t)f Y (r, t). (3) 



By properly ordering the radial channel functions / 7 (t) 
into a column vector F(t) and the evolution operators 
/i 7 and -D 7 . 7 ' into a square matrix H(r, t) we may (in 
the outer region //) rewrite the TDSE of the ejected 
electron of any multi- electron system in the case of single 
ionization as, 

dF 

i—(r,t) = U(r,t)F(r,t) r>b, (4) 

this equation having essentially the form of the one- 
electron TDSE. It is exactly this last equation, no matter 
how the inner region is treated, that allows us to utilize 
any propagation technique in the outer region // of con- 
figuration space, which may have already been applied 
to one-electron ionization. 

On the other hand, in the inner region an eigenstate 
representation of the TD wavefunction will result in a 
TDSE where only two dynamical quantities are needed 
to be provided for the forward propagation in time of the 
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Figure 1: Partition of configuration space for the electron co- 
ordinate. In the inner region / an eigenstate expansion rep- 
resentation of the wavefunction is chosen, while in the outer 
region II a grid representation is considered. 



solution, namely eigenenergies and transition matrix el- 
ements between the system's eigenstates. The key point 
in this case is that the whole information about the ex- 
act nature of the system described in the inner region, 
whether multi-electron or not, is contained in the val- 
ues of the energies and the transition matrix elements 
together with the required selection rules for the tran- 
sitions. Therefore, in a sense, without trying to over- 
simplify, one would expect the matching procedure be- 
tween the two methods (inner region/basis representa- 
tion - outer region/grid representation of the wavefunc- 
tion) to hold regardless of the actual system being multi- 
electron or single-electron in nature. It is for this reason 
we believe the formulation in the present work should 
be readily extendable to complex multi-electron systems. 
The theoretical details and subsequent application will 
be more complicated, due to the multiplicity of ionizing 
channels for the ejected electron in such cases. In the 
following sections we will develop our approach for the 
one-electron atom case in detail thereby laying bare the 
basic concepts of what we believe a novel combination of 
basis set and finite-difference methods. 



III. THE THEORETICAL FRAMEWORK 

In this section we develop a theory for solving the 
TDSE using basis and grid representations in the inner 
and outer regions, respectively. The artificial R-matrix 
division of configuration space into two regions causes 
time-dependent boundary terms to appear in the cor- 
responding TDSEs (in regions I and //) which exactly 
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account for the amount of probability current passing 
through the boundary during the interaction with the 
external field as well as after its turn off. Since the time- 
dependent wavefunction consists of two parts, a careful 
analysis is necessary in order to obtain the physical ob- 
servables of interest such as bound and ionization proba- 
bilities as well as energy and angular information on the 
ejected electron. 

In Sub-Sec. IIII Al we present the calculation of the R- 
matrix eigenstates defined in region / and derive the 
time-evolution equations for a wavefunction expanded 
over the R-matrix eigenstates of the field-free Hamilto- 
nian. In Sub-Sec. IIII Bl we derive the finite-difference 
TDSE governing the radial motion of the ejected electron. 
In Sub-Sec. IIII CI we summarize the calculational pro- 
cedure for the forward in time propagation of the wave- 
function. Finally, in Sub-Sec. IIII Dl we give the formal 
expressions for the calculation of experimental observ- 
ables adapted to our methodology. 

Before proceeding further we first define the inner and 
outer region as shown in Fig. [TJ In region / (defined 
as [0,6]) the TD wavefunction ipj is expanded over the 
eigenstates of the Hamiltonian matrix representation in 
the interval [0,6]. In region II (defined as the interval 
[b, R}) the TD wavefunction is represented by its values 
[ipn(ji, t)] at equidistant grid points r(i) = ih,i = ib, ib + 
1,...,N. 



where the Bloch surface terms Lh and Ld are set out in 
Eqs. (|B2|) and (|B5|) respectively. In these circumstances 
the TDSE over region I is written: 



•Dj(t)]Vi(r,t)- 



+ L d {t)\ i>(r,t), 
(7) 

with < r < 6. This equation is a key one to the method. 
The second term on the right hand side compensates for 
the Bloch terms introduced to make Hi and Dj Hermi- 
tian. Note that it makes a contribution only at r = 6 
and brings into play there ip(r, t) [Eq. (|A3j) ]. a wavefunc- 
tion form which we have defined throughout both regions. 
This term is central to any time propagation scheme in re- 
gion / because it connects the wavefunction form tpi(r, t) 
specific to that region (which may be multi-electronic in 
a more general formulation) with a wavefunction form 
that at r = b represents a single electron and which in 
calculations is obtained from region II. We obtain from 
Eq.0 the evolution equations for the coefficients Cki(t) 
by projection over the states (P k i(r)/r)Yio(8, <fi): 

i-^-Cki(t) = [eui'Skk'Sw + Dki,k'i'(t)] Ck'i'jt) 



- -P u {b)F{{b,t). 



A. R-matrix basis-set TDSE in the inner region 

In the inner region / we define the radial channel func- 
tions fi(r) which are expanded over the R-matrix basis 
set Pki(r) ( defined in appendix IB]) as: 

K 

fi(r,t) = J2 c ki(t)Pki(r), < r < 6, (5) 

k=l 

with the bar at the top indicating that the channel 
function has been obtained by summing over the radial 
Hamiltonian eigenstates of the inner region. Note that 
we have ignored the dependence on the magnetic and 
spin quantum numbers. From the above definition and 
the TD wavefunction [Eq. (|A3j) ] we obtain in the inner 
region /: 

k l p M 

Mr,t)^^C H (f)^ (^), 0<r<6. 

k=l 1=0 

(6) 

The time evolution of the TD wavefunction is now en- 
tirely contained in the coefficients Cki{t). The time evo- 
lution of the Cki{t) is determined by the TDSE. How- 
ever in writing the TDSE we must take care that the 
Hamiltonian and dipole operators which act on ipi(r,t) 
are Hermitian over region I (where ipi(r,t) is only de- 
fined). The Hermitian inner- region Hamiltonian is given 
by Hi = Ha+Lh and the dipole operator by Di = D+Ld, 



The quantity F[(b,t) is defined as: 

i,/ (M )=*M-^ £ K u , Mb ,t), (8) 
ar c ' 

i'=i±i 

where A(t) is the time-dependent field potential in 
the Coulomb gauge (see appendix and Kiy is 
an angular factor given in Eq. (|A5c[) . If the 
coefficient vector C(t) is structured as C T (t) = 
[Cro(t), Cko, Cn(i), ...Cki, Cix(t), CKh{t)\ the 
inner-region TDSE in matrix notation is as follows: 

Cki(t) - -i[H • C] w (t) + l -w k iF{(b,t). (9) 

The amplitudes Wki have been defined as Wki = Pki(b) 
in appendix [Bj The matrix H has the block-triangular 
form of Eq. (|A8|) with the block-diagonal matrices hi 
and the lower- and-upper block matrices Dui (t) having 
matrix elements as: 

(kl\hi\k'l)=e k i6 kk ,, (10) 

(kl\D u±1 \k'l ± 1) = -i^-K n±1 t k wi±i(r), (11) 
c 

where t k i k n±i(r) are matrix elements defined in Eq. 
(IB61) . 
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B. Finite-difference TDSE in the outer region 

In the external region II a grid representation of the 
TD wavcfuction is adopted: 



1=0 



MM) 

r(i) 



Yio{f), b<r(i), (12) 



with i = if,,..,/. The time-dependence of the wave- 
function is represented by the values of the radial 
channel functions on a equidistant discretized grid, 
fl(i,t) = fi(r(i),t) with h = r(i + 1) - r(i),i = i b , 
The grid is defined such that r(i b ) — b and r(I) = R. 
Furthermore by constructing the vector F(t) from the 
values of the radial channels fi(i,t) at the grid points, 
we obtain a vector of length L x I structured as F T (t) = 

[fo(i b ,t),...Jo(I,t)Ji(i b ,t),...J 1 {I,t),...,f L (i b ,t),...f L ^,t 
The FD representation of the TDSE takes the form: 



f l (i,t)=-i[U-F] l (i,t). 



(13) 



In the FD representation of the time-dependent Hamilto- 
nian H(i, t) the entries hi and Du±i(t) are square matri- 
ces of order I— The explicit form of these operators 
depends on the approximation chosen for the derivatives. 
In the present case, the first and the second derivative of 
a function <p(r) are approximated with a 5-point central 
difference scheme as follows: 



2 

E 

j=-2 



hi 



1,2, 



(14) 



with dj chosen so that polynomials of order 4 are dif- 
ferentiated exactly. Given the above, the finite-difference 
approximation of the diagonal operator in the FD Hamil- 
tonian is, 



2 



3=-2 



1(1 + 1) 

2[r(i)] 2 



V(i) 



of boundary conditions imposed on the wavefunction at 
the far boundary. 

Thus some further consideration of the differential op- 
erators involved in the FD representation of the TDSE is 
necessary and we shall shortly see that non-zero function 
values at an inner boundary r — r(ib) — b bring about 
contributions from functions values at points below the 
inner boundary point to the propagation. 

We begin by appreciating that since the FD method 
is a local method, the evaluation of function derivatives 
at any point relics on function values at neighbouring 
points, and which of these come into play depends on 
the approximation chosen for the derivatives, as men- 
tioned earlier. In the FD method the operators are also 
discretized in a similar way to the functions, i.e. as 
0(r, t) = 0(i, t). The action of a non-derivative operator 
n a function is trivial, since 0(r)<f>(r) = 0(i)f(i) at the 
th grid point, but the same is no longer true when oper- 
ators contain derivatives. Then the rule of differentiation 
should be given. The central characteristic of the differ- 
ential operators in the FD method is that values of the 
wavefunctions at neighbouring points are involved in the 
calculation of the derivative function. It is then obvious 
that since the diagonal operators in the finite-difference 
TDSE [Eq. (fl~5|) ] involve the second-order differential op- 
erator (due to the kinetic term) the complete determina- 
tion of the H-F requires knowledge of the /;(«, t) at points 
i = if, — l,ib — 2 since these enter the determination of 
second-order derivatives at points ib and i b + 1 according 
to Eq. (fT4"]) . If the propagation is done in the velocity 
gauge a similar conclusion is reached by considering the 
non-diagonal operators Dui(t). The modified form of the 
TDSE corresponding to a non-vanishing solution on the 
inner boundary is then 

Mht) = -t[H.F],(M) (16) 
+ 5u b [Boi(i b -l,t) + B ol (i b -2,t)} 
+ S iib+1 Bu(i b - l,f), 



■ where 



(15) 

The velocity form of the non-diagonal operator is given 
by: 



,(2) 



Bu(i b -i,t) = - 2^/i(*6 - M) + -ffai^b - i, t) 



(17a) 



_ -iA(t) 



£j=- 3 ir/!±i(i + J\*) 



B Q i{i b - l,i), = - — fi(i b - l,t) + —r-9l(ib - 1,*) 



(17b) 



K, 



«±i- 



B 0l {i b -2,t), 



2h 2 



d 



(i) 



fl(i b -2,t) + — -ffi(* 6 — 2,t) 



The FD form of the TDSE in Eq. |13J) is sufficient 
to propagate the wavefunction in time provided it van- 
ishes at both ends of the spatial grid at all times. This is 
certainly the case when the FD grid has its innermost 
point at the origin. In contrast, in the present case, 
vanishing boundary conditions occur only at the far end 
of the grid (r = R). More specifically, we assume that 
fi(I — l,t) = fi(I, t) = for all I and this forms the set 



(17c) 



and gi (r, i) are given by, 
A(t) 

gi(i,t) = -i + Km +1 )fi+i(i,t)] . 



The elements Kw are given by Eq. (|A5c[) but when I = L 
the term with Km+x) is missing and when I = the 
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term with Kn_x\i IS a l so missing from the corresponding 
equations. The bar on the fi,gi emphasizes that these 
radial function values have been evaluated by use of the 
R-matrix basis set expansion form of the wavefunction in 
region /. 

Eq. (fl6|) is the second (and last!) key equation of the 
method. It does for region II what Eq. (J7J) above did for 
region II. The communication with the solution in region 
/ is provided through the terms involving radial function 
evaluations at two FD points in region / immediately in- 
side the boundary with region II. Although our detailed 
exposition above has centred around one-electron wave- 
functions throughout both regions, it is clear how the 
concept embodied in Eq. (flu) can be extended to handle 
a region / that is multi-electron in character. The crucial 
requirement of such a multi-electron inner region is that 
it must collapse to one-electron character within a few 
FD points of its outer boundary at r = b. Since in multi- 
electron R-matrix calculations anyway the inner region 
must be one-electron in nature by r = b, our additional 
requirement provides no great extra overhead. 



C. Calculational procedure 

Having set out the form of the TDSE in the two regions 
/ [Eq. ©] and II [Eq. ([To]) ] we now present briefly the 
computational procedure involved in the propagation of 
the wavefunction ip(r, t) through one time-step from time 
t to time t + r). 

a. Outer region: calculation of ipu (r, t+r): Assum- 
ing at time t the wavefunction is known throughout the 
inner and outer regions I, II wc first consider the outer 
region II TDSE [Eq. fir?])] . Although there is a wide va- 
riety of methods in the literature we have chosen to em- 
ploy the standard Taylor propagator as prescribed in Eq. 
(IC2[) . The evaluation of the Taylor series terms requires 

the quantities B { °\i b -l,t), B$(i b - 1, t),B$(i b -2,t) 
which bring into play values of the partial waves fi(i — 
2,t),f(i — l,t) evaluated in the internal region at time 
t [Eq. (fTT|) ]. These inner-region partial wave values are 
formed using Eq. 

b. Inner region: calculation of ipi(r,t-\-r): In a sim- 
ilar way as done for the outer region, the propagation of 
the coefficients Cki (t) from time t through one time-step 
to gain their values Cki(t + r) at time t + r is now based 
on the inner-region TDSE in the form of Eq. © and the 
Taylor expansion Eq. (|C2[) . For this evaluation knowl- 
edge of the quantity F{(b,t),l = 0, 1,..,L at time t is 
required. The latter quantity includes the outer-region 
partial wave fi(b,t) and its derivative f[{b,t) evaluated 
on the boundary r = b. Having calculated the coefficients 
Cm (t + t) we can immediately form the wave function 
ipl(r, t + r) according to Eq. ©. 

By this stage the wavefunction is known at time t + r 
throughout regions I and II and we can proceed further 
in time by repeating the above procedure for successive 
time-steps t. 



D. Observables within the dual representation 

In this section we develop the necessary formulation 
for the calculation of observables given the different rep- 
resentation used of the time-dependent wavefunction in 
the inner and outer region (regions / and // respectively). 
These representations are given by Eq. © and Eq. (fT2"j) . 
respectively. Any spatially dependent observable repre- 
sented by the operator 0(r,t) is calculated through the 
standard formula, 0(t) — (ip(r,t)\0(r,t)\ip(r,t)) which 
in our case separates into two pieces. To link with the 
standard experimental setups we assume that any cal- 
culation of the observables is performed for times where 
the external field has vanished. In the following formu- 
las, taking the pulse duration as T, we assume the pro- 
jection time t p such that t p > T. To obtain the pop- 
ulation W n i(tp) in an eigenstate of the physical system 
<Ani(r) = {F n i(r)/r)Yio(f) at time t p , we use the projec- 
tion operator P n i = \<j) n i) {<j>ni\ with the result: 



W nl (t p ) = K^nil/Oz + OF^I/O 



II 



(18) 



with fi(r,t) given by Eq. ((5]) and (a\b)i, (a\b)u denoting 
radial integrations over the inner and outer regions, re- 
spectively. Complete information about the final state 
(ignoring spin variables) is possible by recalling the par- 
tial wave expansion of a continuum electron with asymp- 
totic momentum k = (k, 6k, 4>k), namely: 



(r) =Y,ai mi {k)\F u {r)Y^ mi {k)Y lmi {r), 

Imi 



(19) 



where k = (6 k, 4>k) defines the direction of the photoelec- 
tron with respect to the polarization axis (quantization 
axis), Fki(r) is normalized on the energy scale and the 
amplitudes ai mi (k) are chosen so that the wavefunction 

7Ak _) (r) fulfils incoming spherical wave boundary condi- 
tions. In the present case, where the ionizing target is 
hydrogen and mi — (in the following again we drop 
the mi dependence) we have ai(k) — i l e - lcr d k ) with ai(k) 
the long-range Coulomb phase shift analytically known 
[35l | . Therefore the desired angular distribution is ob- 
tained through the projection operator i\ = ^)(<^k ^| 
which gives: 

dW(e k ,k,t p ) 
dk 



J2 [i F kl\fi)i + (Fki\fi)n] ai(k)Yi (k) 



with dk = k 2 dkdflk the volume element in momentum 
space. Integration of the above formula over the kinetic 
energies (e^ = k 2 /2) results in the photoelectron an- 
gular distribution (PAD), 



dW(e k ,t p ) 
dflh 



dkk 



2 d W{e k ,k,t p ) 
dk 



(20) 



while integration over the photoelectron ejection angles 
(Sk^k) provides the angle-integrated photoelectron en- 



7 



ergy distribution (PES), 
dW(e k ,t p ) 



de k 



{Fki\fi)ii\ 



k=j2^ 



(21) 



Finally, further integration over the photoelectron kinetic 
energies of the last equation results in the total ionization 
probability (yield) at time t p as: 



dW{e k ,t) 
de k 



(22) 



At this point we have completed the present theoretical 
formulation leading to the calculation of the most impor- 
tant experimental observables following the interaction of 
an electromagnetic field with a one-electron atomic tar- 
get in the dipole approximation. 



IV. ILLUSTRATIVE APPLICATION TO 
HYDROGEN 

In the present section we apply our approach to the 
case of ionization of the hydrogen atom by a strong EM 
field. The reasons we have chosen hydrogen are as fol- 
lows: (a) it represents the simplest among the atomic 
systems having just one electron participating in the ion- 
ization process, thus being free from complications that 
may arise from inter-electronic effects in the case of multi- 
electron systems, (b) angular momentum considerations 
are reduced to the minimum level where a simple partial 
wave expansion is adequate to represent the TD wave- 
function throughout the electron's configuration space 
and (c) last but not least ver y re liable methods treating 
one-electron systems [H, H3, 38] are at our disposal for 
a systematic study of the reliability and accuracy issues 
surrounding the present method. In the present appli- 
cation we have chosen an explicit type time-propagator 
based on a Taylor expansion [Eq. (|C1[) ] . In all the calcu- 
lations the order of the propagator was P = 12 and the 
time step r = 1.5625 x 10 -4 a.u. . 



A. Initial state calculation 

We start by calculating the P%o(r), < r < b radial 
function by numerically solving the radial SE for I = 
[Eq. (|B3[) ] within the inner region. The initial state, 
made up of ipi(r, t = 0) and ipn(r, t — 0) in the inner and 
outer region, respectively, is then calculated by means of 
an imaginary time propagation of the field-free versions 
of Eqs. ([9]) and (fT6|) with initial conditions: 

Ckl{t = 0) = 5kl;W, 

fi(i,t = 0) = S l0 S iib P 10 {b). 

It is important to emphasize here, that the R-matrix 
eigenstates do not actually represent the eigenstates of 
the system, instead they only serve as a complete basis 




Figure 2: (Color online) Initial state wavefunction calculated 
by an imaginary field-free time propagation. Fhe spatial step 
was h = 0.29 a.u. and the time step was r = 1.5625 x 10~ 4 
a.u. The order of the Taylor propagator was 12. Black curve 
represents the trial wavefunction while the red (solid) and the 
green (dashed) curves represent the converged wavefunction 
for the inner and outer regions, respectively. 



for the representation of the physical state exclusively in 
the interval [0, R]. The B-splines basis used consisted 
of nf, — 57 basis functions of order — 9. The knot- 
sequence is chosen to be linear with discretization step 
equal to that of the outer region spatial step h = 0.29 
a.u.. In Fig. [2]we plot the squared amplitude of the radial 
part of the R-matrix basis for n = 1, 1 = 0, Pio(r) (black 
curve) and the state Pi s (r,t — > oo) = r(Yio\ip(r, t —> 
ioo)), < r < R as converged after the imaginary-time 
field-free propagation. Black-solid and red-dashed curves 
represent the inner (/) and outer (II) region values. The 
inner boundary has been set at b — 14.5 a.u. while the 
outer boundary at R = 174 a.u. The R-matrix eigenstate 
Pio(r) has zero-derivative on the boundary (due to the 
chosen boundary condition P' kl (b) — 0) while the imagi- 
nary time propagation has converged to the state Pi s (r) 
with non-vanishing derivative on the boundary as actu- 
ally is the case for the ground state of hydrogen. Our 
initial state has the following form: 



^(r,0) = -Y 10 (f) 
r 



J2 k Cfeo(0)Pfeo(r), region / 



/oM), 



region // 



(23) 



B. Real time propagation 

Having obtained an accurate initial state [Eq. (f2"3"|) ]. 
through imaginary time propagation, we proceed to the 
propagation of the TDSE in the outer/inner region in 
the presence of an external EM field. The EM field cho- 
sen was linearly polarized along the z-axis with vector 
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Figure 3: (Color online) Hydrogen bound state population 
within the (0, 14.5) a.u. region after irradiation by an external 
EM field (see text for the field details). Curves represent the 
present mixed method (BS/FD) as well as standard finite- 
difference (FD) and eigenstate expansion (BS) methods. 



Figure 4: (Color online) Absolute square of the hydrogen par- 
tial wave I = 1 (\fi{r,t p )\ 2 ) at t p = 157.1 a.u. after irradiation 
by an external EM field (see text for the field details). Curves 
represent results with the present mixed method (BS/FD) as 
well as with the standard finite-difference (FD) method. 



potential: 

A(t) = A sin 2 (^t) smut, (24) 

where to = 2t:/Tq is the field frequency and T = uTq 
the pulse duration (n being the number of cycles con- 
tained in the pulse and To the field period). To test our 
approach we calculated the bound state population using 
three different methods. First, within the present method 
(BS/FD), we obtained the population in the [0, b] region 
by simply summing (at a sufficiently long time t p ) over 
all R-matrix eigenstates as: 

#(*p)=£|Cw(tp)| a . (25) 

kl 

The second method consisted of the standard finite- 
difference (FD) approach over both regions I and 11 i.e. 
over the whole range [0, R] (with the same spatial and 
time step) thereby invoking no division of the electron's 
configuration space. Formally, within the present method 
this is equivalent to setting b — 0. We calculated the ion- 
ization probability P c as: 

r(i)>6 

Pc(t P )= l^,i P )| 2 , (26) 
a 

where t p was chosen sufficiently large so that all the out- 
going components of the ionized wavepacket were able 
to travel beyond the chosen distance b. With no ab- 
sorbing potential present we always have for the bound 
state probability Pb(t) = 1 — P c (t). When an absorbing 
potential is present then the bound state probability is 
obtained as: 

r(i)<b 

A(tp)= E I/<(Mp)| 2 - (27) 

a 



Finally, we performed calculations using a standard ba- 
sis set (BS) to span the whole range [0, R] with again no 
division of configuration space. This is formally equiv- 
alent to setting b = R. We obtained the bound state 
population by summing only over the bound part of the 
spectrum: 

n CBS) (* P )=£|CW(* P )| 2 , e kl <0. (28) 

kl 

In Fig. [3] we show the bound state population of hydro- 
gen as a function of time when the atom is irradiated by 
a pulse of central frequency w = 0.8 a.u. (21.769 eV), to- 
tal duration of 10 cycles (T = 10T = 10 x 7.854 = 78.54 
a.u.) and peak intensity Iq = 10 14 W/cm 2 . For the 
BS calculation the bound state population is calculated 
at the end of the pulse (10 cycles) and no further field- 
free propagation of the wavefunction is required since the 
population distribution remains unchanged. In the case 
of the FD and BS/FD calculations the propagation is 
extended for a further 10 cycles (field- free propagation) 
after the end of the pulse until a sufficiently large part 
of the wavepacket has passed the artificial boundary at 
r = b = 14.5 a.u.. Given the photon frequency, the hy- 
drogen ionization potential and the rather modest field 
intensity, we expect the dominant partial wave in the 
outer region to be the I = 1 partial wave with the elec- 
tron's kinetic energy peaked around e& ~ 0.8 — 0.5 = 0.3 
a.u. (8.16 eV). By assuming an outgoing wavepacket with 
central energy of 0.3 a.u. (thus of velocity k = 0.7746 
a.u.) we can estimate that 10 cycles of field-free propaga- 
tion is sufficient for our purposes. In connection with this 
latter point note that this wavepacket travels a distance 
of 14.5 a.u. in approximately 2.5 field-cycles. This is 
why the FD and BS /FD bound state populations exhibit 
a time-delay compared to the BS bound state population. 
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The maximum angular momenta allowed was L — 3. Re- 
sults remained practically unchanged against further in- 
crease in angular momentum. We have performed similar 
calculations with peak intensities Iq = 10 15 W/cm 2 and 
found similar results with analogous agreement between 
the BS/FD, FD and BS bound state populations. 

In Fig. |4] the values of \fi(r, t p )\ 2 , 1 — 1 are plotted 
as calculated with the present BS/FD and the standard 
FD method at t p = 20T = 157.1 a.u.. In region I 
(within 14.5 a.u. of the nucleus) the partial wave func- 
tion (|/i(r, t p )\ 2 ) was obtained from Eq. ([5]). In region 77 
(from 14.5 a.u out to 174 a.u.) the values of |/i(r, t p )\ 2 
come directly from the propagation of the outer-region 
TDSE [Eq. (fj"6])]. Similarly for the FD calculation we 
obtained |/i(r, t p )\ 2 by solving Eq. (fl~3"[) over the whole 
range [0,R]. The figure displays excellent agreement be- 
tween such results from the present (BS/FD) method and 
the standard FD method. We have chosen to plot only 
the I = 1 partial wave since this is the dominant out- 
going channel with all other partial wave channels being 
an order of magnitude lower. This observation simplifies 
the analysis of the physics involved in the process. We 
briefly elaborate on this plot. The peak probability for 
the travelling wavepacket appears around ~ 92 a.u. with 
a much smaller secondary peak inside region /. In an en- 
ergy representation of the wavepacket, the large peak is 
associated with the continuum states contribution while 
the second peak is related to the bound states contri- 
bution. Whereas the bound contribution is trapped in 
the inner region the outgoing component (corresponding 
to the continuum spectrum) travels a distance of about 
r ~ v x 15T = 0.7746 x 15T ~ 91 a.u. which is rather 
close to the maximum of the wavepacket probability in 
the plot. We have allowed 15 cycles of travelling time 
for the wavepacket since significant ionization only takes 
place around the maximum of the applied pulse which 
occurs at approximately 5 cycles after the turn-on. 

In Fig. [5] the bound state population of hydrogen is 
shown after exposure to an EM field of central frequency 
lo = 0.35 a.u. (9.524 eV), total duration 10 cycles and 
peak intensity Iq = 10 14 W/cm 2 . Since the photon en- 
ergy is comparable to the energy gap (~ 10.277 eV) be- 
tween the ground and first excited states (2s, 2p) an ap- 
preciable population in these excited states appears at 
the end of the pulse. At the end of the pulse we ob- 
tained from the BS calculation a value P g — 0.7368 for 
the ground state probability ; a value P e = 0.2104 for the 
total population in all the excited states (ejy < ) and 

thus a total bound state probability of P b = 0.947318 
[Eq. ([2"8")) ]. The bound state probability as a function of 
time is shown in the figure (blue line) . We have also per- 
formed a FD calculation (with no absorbing boundary 
present) and calculated the probability within the region 
[29, R] a.u. using Eq. ^ and P b = 1 - P c . We chose 
a box with R = 522 a.u. to prevent reflection of the 
wavepacket at the outer boundary over the time interval 
of interest. In this case the calculated bound state prob- 
ability for the FD method is given by the green curve 




m—m BS/FD, Eq. (28) 
•— • BS/FD, Eq. (25) 
0- -e FD, Eq. (27) 
BS. Eq. (28) 



Figure 5: Hydrogen bound state population within the (0, 29) 
a.u. region after irradiation by an external EM field of pho- 
ton frequency of uo — 0.35 a.u. and peak intensity I — 
10 14 W/cm 2 . Curves represent the present mixed method 
(BS/FD) as well as standard finite-difference (FD) and eigen- 
state expansion (BS) methods. 



(empty cycles). Next, we applied the present method 
(BS/FD) for b = 29 a.u. and R = 552 a.u. To maintain 
the same accuracy in the calculations in the inner region 
we increased the number of B-splines basis members to 
rib = 108. To compare with the BS and FD calculations 
we obtained the various probabilities as follows: Black 
curve (filled squares) in the figure was calculated using 
Eq. (|25|) which includes a summation only over those 
R-matrix eigenstates that have negative energies such 
that €k < [equivalent to Eq. |2"5[)]. This curve follows 
closely the bound state probability calculated using the 
BS method. If in Eq. ([25]) we include all R-matrix states 
then the probability enclosed in region / is given by the 
red curve (filled circles) and matches perfectly with the 
bound state probability from the FD calculation. Similar 
evaluation through Eq. ([26]) and Pf, = 1 — P c with the 
BS/FD method results in practically the same curve and 
verifies the equality of results obtained using Eq. (j2"5|) 
and Eq. ([26]) . In other words any increase/decrease of 
probability within region I is matched by an equal de- 
crease/increase of probability within region II. 

Finally in Fig. [5] we have calculated the photoclec- 
tron energy spectrum up to about 20 eV kinetic energy 
of the ejected electron. In the hydrogen case, although 
the analytical solutions for the bound and the continuous 
spectrum are available, the numerical calculation of the 
eigenstates proves more advantageous for the evaluation 
of the necessary integrals. In expression (|21|) we may 
use an asymptotic expansion for the Coulomb functions 
Fu (r) [23, [40] provided that (a) the evaluation is per- 
formed at times where the outgoing part of the electron 
wave packet has travelled sufficiently far away from the 
residual system (b) the projection operator is constructed 
either from Coulomb wavefunctions or plane waves de- 
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Figure 6: Hydrogen photo-electron energy spectrum after ir- 
radiation by an external EM field of photon frequency of 
w = 0.35 a.u. (9.52 eV) and peak intensity I = 10 14 W/cm 2 . 
Curves represent the present mixed method (BS/FD) and the 
eigenstate expansion (BS) methods. 

pending on the chosen projection time (t p ) and (c) the 
inner-region contribution is ignored since it is only the 
bound part of the wavepacket that still remains there as 
time grows. The results can easily be checked by tracing 
their convergence in time. A detailed discussion of this 
approach, very well suited to our approach, can be found 
in ref. [39|]. The solid black curve represents the result 
of a BS calculation while the dashed red curve the result 
of the present calculation. Had we used a larger box and 
a finer mesh for the outer region we would be able to 
calculate even higher (in energy) the corresponding PES 
for a full comparison with the BS calculation, but this is 
not the purpose of the present work. 



V. CONCLUSION AND PERSPECTIVES 

In conclusion a new ab initio time-dependent method 
for the treatment of the single-electron ionization of 
atomic and molecular systems under an external electro- 
magnetic field has been set out. It has been developed 
in detail for systems that are single-electron throughout 
and applied to the simplest case, namely the hydrogen 
atom. The method is based on the division of the con- 
figuration space of the ejected electron into two regions 

I and II. In region / (which may be multi-electronic) 
the time-dependent wavefunction is expanded on the ba- 
sis of R-matrix eigenstates and propagated through the 
time evolution of the expansion coefficients. In region 

II a grid-representation of the time-dependent wavefunc- 
tion is adopted and a finite-difference technique is em- 
ployed for the representation of the operators. In both 
regions the chosen time-propagator in illustrative calcu- 
lations is a high-order explicit Taylor propagator. The 
key point in the present method is the time-dependent 



matching conditions that the inner (region /) and outer 
(region //) wavefunctions should simultaneously satisfy 
at each time step. Although these matching conditions 
have been developed here for an explicit time-propagator, 
the methodology can also be applied for implicit time 
propagators. The present work represents an impor- 
tant step towards the implementation of such a method- 
ology in multi-electron systems (atomic and molecular) 
where the full advantage of the R-matrix technique can 
be taken into account. The straightforward extension of 
the present approach to the case of a truly multi-electron 
system is discussed in Sec. Ulland is currently the subject 
of our work. In addition to our fundamental interest in 
gaining an ab-initio description of multi-electron systems 
under strong laser fields the present work is mainly moti- 
vated by the development of sources of short wavelength 
laser light residing well into the VUV or soft XUV regime 
(HOHG/FEL sources). In contrast to long wavelength 
laser light, the light from such sources tends to inter- 
act directly with more than just a single electron and is 
able to probe directly the innermost electrons of multi- 
electrons systems thus making the development of new 
suitable theoretical methods a necessary and formidable 
task. 

The authors gratefully acknowledge discussions with 
Dr. Michael Lysaght, Dr. Hugo Van der Hart and 
Prof. P. G. Burke. The present work has been supported 
through a Marie Curie Intra-European Fellowship under 
contract TDRMX-040766 awarded to LAAN and by the 
UK Engineering and Physical Sciences Research Council. 

Appendix A: TDSE OF 
SINGLE-ACTIVE-ELECTRON ATOMIC 
SYSTEMS OVER A SPHERICAL HARMONIC 
BASIS 

The field-free SAE Hamiltonian Hq reads, 

H = -±V 2 + V(r), (Al) 

with the potential V(r) equal to —Z/r for a purely hy- 
drogenic system (of Z atomic number). Alternatively 
V(r) could be constructed as a model or Hartree-Fock 
potential. The TDSE of the system in an external time- 
dependent radiation field, E(f) is written as: 

ij^(r> t) = [Ho(r, t) + D(v, t)] ^(r, t), (A2) 

with ip(r,t) the system wavefunction and D(r,t) the in- 
teraction operator between the system and the external 
field, in the dipole approximation. In our present nu- 
merical implementation we choose a spherical coordinate 
system for the active electron. We represent the angular 
variables in a basis of spherical harmonics and write the 
wave function as, 

^M) = £ E ^^lWM), (as) 

1=0 m=-l 



11 



where the spin- variables of the wavefunctions are ignored. 
In an actual calculation we must truncate the spherical 
harmonics expansion at some maximum value L. In the 
remaining formulas we abbreviate the truncated double 
summation by J2i m - 

The time-propagation of the wave function proceeds in 
spherical coordinates as follows. Substituting Eq. (|A3[) 
in Eq. (IA2j) and projecting onto the spherical harmonic 
basis Yi m (f) we obtain the following coupled differential 
equations for the radial channel functions as, 

d - 

ig^fim(r,t) = hi(r)fi mi (r,t)+22 Am,,i'ni;(' , |i)/i'm;('' ) t). 

I'm', 

(A4) 

For the special case of linearly polarized light along 
the z-axis and in dipole approximation the radial time- 
evolution operators are given by 



hi(r) 



1 d 2 + 



+ V{r), (A5a) 



Di mi ;i'm'(r,t) = -i S mu m,,Kiu(mi)tw(r), (A5b) 

i c 



I 2 - m 2 

Ku>(mi) = S u±1 J > 2 _ J , 



(A5c) 



with /> = max(l,l') and ty> the radial dipole operator. 
The time-dependent radial dipole operator is given as, 



in>(r) 



d_ 

dr 



(I -I') 



(A6) 



in the velocity form. The quantity A(t) = — c J Q dt'E(t') 
represents the held potential in the Coulomb gauge. 
Within the present context the interaction operator cou- 
ples atomic states of equal magnetic quantum number, 
hence we drop the dependence on mi in the subsequent 
formulation. 

By properly arranging the radial channel functions /; 
according to their angular momentum label we form the 
radial vector wavefunction F. In this case the matrix 
representation of the TDSE [Eq. (| A4[) ] is written as, 



F(t) = -»H(t)F(t), 



(A7) 



where F = dF(t)/dt and 



H(r,t) 



ho Doi 
D w h 2 D12 

^21 h 







hi,-x Dl-i,l 

£>L,L-1 h L 



Appendix B: R-MATRIX EIGENSTATES IN THE 
INNER REGION 



1. Hamiltonian operator in the inner region 

In the inner region [0, b] the radial wavefunctions 
fi (r, t) are expanded over the eigenstates of the radial 
Hamiltonian: 

h = hi+L h l = 0,l,..,L. (Bl) 
with Lh the radial Bloch operator, 

l h = \s{r-b)± (B2) 

and hi given by Eq. (|A1[) . The eigenstates of the R- 
matrix Hamiltonian operator hi are uniquely determined 
if we set the boundary conditions needed to be fullhled 
at the boundaries r = and r = b. In the present case on 
physical considerations we take all solutions to vanish at 
the origin while at r = b the solutions take non- vanishing 
values. This choice makes the radial R-matrix operator 
Hermitian over the inner region [0, b]. Therefore for each 
value of the angular momentum we solve the following 
eigenvalue problem: 



hiP ki {r) =e k iP H (r) 



1 = 0,1, 



(B3) 



where k is an integer labelling the eigenstate. The above 
eigenvalue differential equation is transformed to solv- 
ing a matrix diagonalization problem by employing a B- 
spline basis set of size rib, order fc& [4l| for the represen- 
tation of the solutions Pu (r) in region I: 



Pki(r) 



J=2 



(W) n (fe 6 ) 



< r < b. 



(B4) 



In the expansion the first B-spline [B[ hb \r)] is excluded 
in order to conform to the boundary condition at the 
origin Pjy(O) = 0. Note that by definition of the B- 
splincs the amplitude of the eigenstates on the boundary 
[wki = Pki(b)] is simply the last coefficient in the ex- 
pansion, namely, wm — All required integrals are 
evaluated, with the Gaussian quadrature rule, to machine 
accuracy. 

For each partial wave I = 0, 1, ...,L the solutions con- 
stitute an orthonormal basis with rib — 1 members, 

rib — i 

\ p ki)(Pki\ = i, (Pki\p k >i) = s kk >, 

k=l 

with real eigenvalues e^. 

(A8) 2. Dipole operator in the inner region 

While the velocity form of the radial dipole opera- 
tor includes a first-order derivative term [Eq. (|A6[) ] 
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which taken together with the non-vanishing values 
of the eigenstates Pki(r) at the boundary b makes it 
non-hermitian.We can make this operator Hermitian by 
adding the dipole Bloch operator for the first-order 
derivative in a similar way as done for the field-free 
Hamiltonian hi. Thus if we define the dipole velocity 
operator in region I as: 



L d = -S(r - b)cosd 

we find for the radial velocity operator: 
fb 



kl.k'V 



drP kl (r) 



tu.k'i' ~ 2% - b ) 



(Wo) 



P k n>{r). (B6) 



Appendix C: TAYLOR PROPAGATOR 

The forward evolution of a time-dependent function 
P(t) from a time t to a time t + r by the time-step r, can 
be approximated by the Taylor expansion [34| : 



F(t 



p 

p=0 



'dtp 



(CI) 



with r = t„+i — in, n = 0, 1, N and a p = t p /p\. The 
above propagation scheme consists of an explicit one-step 
scheme of order P. 

When the evolution equation for the F(t) is known 
as F(t) = —iH(t)F(t) the above expression can also be 
obtained as the P— order expansion of the evolution op- 
erator exp(— iH(t)r): 



F(t + r) = e- l X +Tdt ' H (^F{t) 



(C2) 



The above approximate expressions for the time evolu- 
tion assumes that the characteristic time evolution of the 
Hamiltonian H(t) is much larger than the time step r. In 
other words, the Hamiltonian operator is assumed con- 
stant within the interval [t, t+r] evaluated at time t. Fur- 
thermore, in this summation higher-order time deriva- 
tives (H(i), H(t), ...) of the Hamiltonian operator have 
been dropped, a procedure very well justified for the elec- 
tric field strengths used in this work. 
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